数値微分と非線形最小二乗法

In [54]:
%matplotlib inline
import numpy as np
from matplotlib.pyplot import *

pi=np.pi
sin=np.sin
cos=np.cos
sqrt=np.sqrt;
vec = lambda x: x[:,np.newaxis]
array = np.array
linspace = np.linspace

非線形モデルを扱うとモデルを偏微分しなくてはいけないことがある． ただ，モデルが複雑だと偏微分計算がとても面倒くさかったりする． 元のモデルからマシンパワーによって数値的に偏微分を求めることができるらしいのでメモ

1.数値微分

関数$f(x)$の微分の定義は以下の式である．

$$\frac{f(x)}{dx}=\lim_{\Delta x \to 0} \frac{f(x-\Delta x)-f(x)}{\Delta x}$$

数値微分は$\Delta x$にプログラム上で小さい値を代入して定義の通りに計算する方法である．

In [55]:
def fprime(fn):
    def fp(x):
        accr=0.001
        return (fn(x+accr) - fn(x))/accr
    return fp

実際にsinを数値微分してcosを得る．

In [56]:
th=np.linspace(-pi,pi,100)
f=lambda x : sin(x)
sin_prime=fprime(f)

sin_prime_val=[]
for t in th:
    sin_prime_val.append(sin_prime(t))

figure()
plot(th,sin(th),label='sin')
plot(th,cos(th),'r',label='cos')
plot(th,sin_prime_val,'g--',label='numerical dif')
grid()
legend()
Out[56]:
<matplotlib.legend.Legend at 0x1eea0cb9948>

2.非線形最小二乗法

2.1概要

非線形最小二乗法は例えばGPSで用いられる． 具体的には，「利用者の位置」を「衛星の位置」「衛星-利用者間の距離」を用いて求める問題である． 今回は衛星が3機の場合を考える．

In [57]:
#sat pos
sx1=vec(array([2,2]))
sx2=vec(array([3,1]))
sx3=vec(array([0,1.5]))
sx=np.c_[sx1,sx2,sx3]
In [58]:
#obs
r1=2
r2=1.7
r3=2
y=np.c_[r1,r2,r3]
In [59]:
def plot_implicit(fn,Chi=0,xrange=[-1,1],yrange=[-1,1],resol=100,label='__nolabel'):
    """
    plot implicit function fn(x,y)=Chi (Chi=const.)
    """
    x=np.linspace(xrange[0],xrange[1],resol)
    y=np.linspace(yrange[0],yrange[1],resol)
    x,y =np.meshgrid(x,y)
    z=fn(x,y)
    if label =='__nolabel':
        contour(x,y,z-Chi,[0])
    else:
        contour(x,y,z-Chi,[0],label=label)
    
#plot

def drawModel(sx,r):
    return lambda x,y: (sx[0]-x)**2 + (sx[1]-y)**2 -r**2
xrange,yrange=(-2,5),(-1,4)
figure()
plot(sx[0,:],sx[1,:],'ob',label='satellite')
plot_implicit(drawModel(sx1,r1),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx2,r2),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx3,r3),xrange=xrange,yrange=yrange)
grid()
axis('equal')
legend()
Out[59]:
<matplotlib.legend.Legend at 0x1eea0dcd108>

図を見ると，3つの円の交点付近が利用者の位置であるように思われる． その感覚を数学的に考えようというのが非線形最小二乗法である．

2.2ニュートン法

観測値である距離と，利用者の位置の間には理論的にこんな関係があるはずだ，という理論式を考える．この理論式のことを観測モデルと呼ぶ(個人的に呼んでいるだけかもしれない)．

GPSの場合は，利用者の位置を$x,y$(未知数)，i番目の衛星の位置を$x_i,y_i$,i番目の衛星と利用者の距離を$r_i$とし，$r_i$が観測として得られるとすれば，以下の観測モデルが得られる．

$$r_i=\sqrt{(x_i-x)^2 + (y_i-y)^2}$$

また，観測モデルを行列形式で書くと以下のようになる．

$$\begin{eqnarray} \left( \begin{array}{c} r_1\\r_2\\r_3 \end{array} \right)&=& \left( \begin{array}{c} \sqrt{(x_1-x)^2 + (y_1-y)^2} \\ \sqrt{(x_2-x)^2 + (y_2-y)^2} \\ \sqrt{(x_3-x)^2 + (y_3-y)^2} \end{array} \right) \\ \iff \left( \begin{array}{c} r_1\\r_2\\r_3 \end{array} \right)&=& \left( \begin{array}{c} f_1(\boldsymbol{x}) \\f_2(\boldsymbol{x}) \\f_3(\boldsymbol{x}) \end{array} \right)\\ \iff \boldsymbol{y} &=& \boldsymbol{F}(\boldsymbol{x}) \end{eqnarray}$$

変数置き換えは以下のとおり． $$\begin{eqnarray} \boldsymbol{x}&=& \left( \begin{array}{c} x\\y \end{array} \right)\\ f_i(\boldsymbol{x})&=&\sqrt{(x_i-x)^2 + (y_i-y)^2}\\ \boldsymbol{y} &=& \left( \begin{array}{c} r_1\\r_2\\r_3 \end{array} \right)\\ \boldsymbol{F} &=& \left( \begin{array}{c} f_1(\boldsymbol{x}) \\f_2(\boldsymbol{x}) \\f_3(\boldsymbol{x}) \end{array} \right) \end{eqnarray}$$

In [60]:
# model
def obsModel(sx):
    return lambda usrx : np.sqrt((sx-usrx).T@(sx-usrx))[0,0]
f1=obsModel(sx1)
f2=obsModel(sx2)
f3=obsModel(sx3)
F=vec(np.array([f1,f2,f3]))

非線形最小二乗法の概念をイメージ図を用いて説明する． 観測モデル$\boldsymbol{y}=\boldsymbol{F}(\boldsymbol{x})$をイメージするために，横軸を$\boldsymbol{x}$,縦軸を$\boldsymbol{F}(\boldsymbol{x})$にとるような本来描けないグラフが仮に以下のように描けたとする．

In [61]:
#### test
f=lambda x : sqrt(x)+1
finv = lambda y : (y-1)**2
fprime = lambda x :1/(2*sqrt(x))
testx=np.linspace(0,10,100)
testy=f(testx)

figure()
plot(testx,testy,'k',label='y=f(x)')
xlim((0,10))
ylim((0,4))
legend(loc=5)
Out[61]:
<matplotlib.legend.Legend at 0x1eea0e70208>

ここで，観測値$\boldsymbol{y}$が得られたとき，観測モデル$\boldsymbol{F}(\boldsymbol{x})$が観測値$\boldsymbol{y}$を出力するような$\boldsymbol{x}$を知りたい．

そのような$\boldsymbol{x}$を$\boldsymbol{\hat{x}}$とする．

この問題で困っていることは，観測モデル$\boldsymbol{F}(\boldsymbol{x})$が複雑だから観測モデルを$\boldsymbol{x}=\boldsymbol{F}^{-1}(\boldsymbol{y})$の形に変形することが難しく，その結果$\boldsymbol{\hat{x}}$を求めることが難しいということである．

In [62]:
y=3

figure()
plot(testx,testy,'k',label='y=f(x)')
plot([finv(y),finv(y),0],[0,y,y],'om--')
xlim((0,10))
ylim((0,4))
legend(loc=5)
text(0,y,'y=f(xhat):given',va = 'bottom',fontsize=18)
text(finv(y),0,'xhat : difficult to get',va = 'bottom',fontsize=18)
Out[62]:
Text(4, 0, 'xhat : difficult to get')

そこで，ニュートン法では以下のような考え方をする．

  • まず適当に解を仮定する．仮定した解を$\boldsymbol{x}_0$とおく．
  • このとき，知りたい解$\boldsymbol{\hat{x}}$と仮定解$\boldsymbol{x}_0$の間には$\boldsymbol{x}_0 + \Delta \boldsymbol{x} = \boldsymbol{\hat{x}}$の関係がある．
  • もし$\boldsymbol{F}(\boldsymbol{x}),\boldsymbol{F}(\boldsymbol{x}_0),\boldsymbol{x}_0,\boldsymbol{y}$を用いて$\Delta \boldsymbol{x}$を推定することができれば知りたい解$\boldsymbol{\hat{x}}$が得られる．

このとき，鍵になるのは以下の考え方である．

  • $\Delta \boldsymbol{x}$を正確に求められなくても，近似解が求められれば，$\boldsymbol{x}_0 + \Delta \boldsymbol{x}$は知りたい解$\boldsymbol{\hat{x}}$に近づく．
  • $\boldsymbol{x}_0 + \Delta \boldsymbol{x}$を新たな仮定解$\boldsymbol{x}_0$として，$\Delta \boldsymbol{x}$の近似解を求めて・・・を繰り返せばいつかは知りたい解にたどり着く．
In [63]:
x0=1
figure()
plot(testx,testy,'k',label='y=f(x)')
plot([x0,x0,0],[0,f(x0),f(x0)],'ob--')
plot([finv(y),finv(y),0],[0,y,y],'om--')
plot(linspace(x0,finv(y),10),np.ones((10,1)))
xlim((0,10))
ylim((0,4))
legend(loc=5)
text(x0,0,'x0:given',va = 'bottom',fontsize=18)
text(0,f(x0),'f(x0):given',va = 'bottom',fontsize=18)
text(0,y,'y=f(xhat):given',va = 'bottom',fontsize=18)
text(finv(y),0,'xhat=x0+dx',va = 'bottom',fontsize=18)
text(x0+0.1,1,'dx:target',va = 'bottom',fontsize=18)
Out[63]:
Text(1.1, 1, 'dx:target')

$\Delta \boldsymbol{x}$の近似には以下の式を用いる．

$$\begin{eqnarray} \boldsymbol{F}(\boldsymbol{x}_0 + \Delta \boldsymbol{x}) &\fallingdotseq& \left. \frac{\partial \boldsymbol{F}(\boldsymbol{x})}{\partial \boldsymbol{x}} \right|_{\boldsymbol{x}=\boldsymbol{x}_0}\Delta \boldsymbol{x} +\boldsymbol{F}(\boldsymbol{x}_0) \\ \iff \boldsymbol{y} &\fallingdotseq& \left. \frac{\partial \boldsymbol{F}(\boldsymbol{x})}{\partial \boldsymbol{x}} \right|_{\boldsymbol{x}=\boldsymbol{x}_0}\Delta \boldsymbol{x} +\boldsymbol{F}(\boldsymbol{x}_0) \end{eqnarray}$$

左辺については以下の等式から得られる． $$\boldsymbol{y}=\boldsymbol{F}(\boldsymbol{\hat{x}})=\boldsymbol{F}(\boldsymbol{x}_0+\Delta \boldsymbol{x})$$ また，左辺が示す値は図上ではマゼンタの横線に対応する．

右辺について． 偏微分の項は図的なイメージでいくと関数$\boldsymbol{F}(\boldsymbol{x})$の$\boldsymbol{x}=\boldsymbol{x}_0$における接線の傾きを表す． なので，右辺全体が示す値は図上ではシアンの横線に対応する．

図上で右辺と左辺を比較すると，線形近似臭が感じ取れる．

In [64]:
t=np.arange(0,6)
tngline = lambda x : fprime(x0)*(x-x0)+f(x0)
yt= tngline(t)

figure()
plot(testx,testy,'k',label='y=f(x)')
plot(t,yt,'r',label='y=f\'(x0) (x-x0) +f(x0)' )
plot([x0,x0,0],[0,f(x0),f(x0)],'ob--')
plot([finv(y),finv(y),0],[0,y,y],'om--')
plot([0,finv(y)],[tngline(finv(y)),tngline(finv(y))],'oc--')
xlim((0,10))
ylim((0,4))
legend(loc=5)
text(x0,0,'x0',va = 'bottom',fontsize=18)
text(0,f(x0),'f(x0)',va = 'bottom',fontsize=18)
text(0,y,'y=f(x+dx)',va = 'bottom',fontsize=18)
text(finv(y),0,'xhat',va = 'bottom',fontsize=18)
text(0,tngline(finv(y)),'f(x+dx) Approx.',va = 'bottom',fontsize=18)
Out[64]:
Text(0, 3.5, 'f(x+dx) Approx.')

ここで，イメージから実際の計算に戻る． やりたいことは，以下の式を$\Delta \boldsymbol{x}$について解くことである．

$$\begin{eqnarray} \boldsymbol{y} &=& \left. \frac{\partial \boldsymbol{F}(\boldsymbol{x})}{\partial \boldsymbol{x}} \right|_{\boldsymbol{x}=\boldsymbol{x}_0}\Delta \boldsymbol{x} +\boldsymbol{F}(\boldsymbol{x}_0) \end{eqnarray}$$

それにあたって，まず以下の式を考える． $$\frac{\partial \boldsymbol{F}(\boldsymbol{x})}{\partial \boldsymbol{x}}$$ この式はヤコビアンと呼ばれていて，通常$\boldsymbol{J}$で表す．ヤコビアン$\boldsymbol{J}$は具体的には以下の要素を持つ． $$ \boldsymbol{J}= \frac{\partial \boldsymbol{F}(\boldsymbol{x})}{\partial \boldsymbol{x}}= \left( \begin{array}{c} \frac{\partial f_1(\boldsymbol{x})}{\partial{x}} & \frac{\partial f_1(\boldsymbol{x})}{\partial{y}}\\ \frac{\partial f_2(\boldsymbol{x})}{\partial{x}} & \frac{\partial f_2(\boldsymbol{x})}{\partial{y}}\\ \frac{\partial f_3(\boldsymbol{x})}{\partial{x}} & \frac{\partial f_3(\boldsymbol{x})}{\partial{y}}\\ \end{array} \right)$$

ヤコビアンが計算できれば，通常の最小二乗法と同様に擬似逆行列を用いて$\Delta \boldsymbol{x}$は以下で与えられる．

$$ \Delta \boldsymbol{x} = \left( \boldsymbol{J}_0^\mathrm{T}\boldsymbol{J}_0 \right)^{-1}\boldsymbol{J}_0^\mathrm{T} (\boldsymbol{y} - \boldsymbol{F}(\boldsymbol{x}_0)) $$

ただし，$\boldsymbol{J}_0=\left. \boldsymbol{J} \right|_{\boldsymbol{x}=\boldsymbol{x}_0}$である．

2.2 解析的にヤコビアンを求めて解く

ここからがこのノートの本題で，結局のところヤコビアンをどうやって求めるか？という話である．

一つの方法として，真っ向から偏微分を計算するという選択がある．

前述の通りGPSの観測モデルは以下で与えられる． $$\begin{eqnarray} \boldsymbol{F} &=& \left( \begin{array}{c} f_1(\boldsymbol{x}) \\f_2(\boldsymbol{x}) \\f_3(\boldsymbol{x}) \end{array} \right)\\ \boldsymbol{x}&=& \left( \begin{array}{c} x\\y \end{array} \right)\\ f_i(\boldsymbol{x})&=&\sqrt{(x_i-x)^2 + (y_i-y)^2}\\ \end{eqnarray}$$

この観測モデルに対してヤコビアンは以下で与えられる． $$\boldsymbol{J}=\left( \begin{array}{c} \frac{-(x_1-x)}{\sqrt{(x_1-x)^2+(y_1-y)^2}} &\frac{-(y_1-y)}{\sqrt{(x_1-x)^2+(y_1-y)^2}}\\ \frac{-(x_2-x)}{\sqrt{(x_2-x)^2+(y_2-y)^2}} &\frac{-(y_2-y)}{\sqrt{(x_2-x)^2+(y_2-y)^2}}\\ \frac{-(x_3-x)}{\sqrt{(x_3-x)^2+(y_3-y)^2}} &\frac{-(y_3-y)}{\sqrt{(x_3-x)^2+(y_3-y)^2}}\\ \end{array} \right)$$

この方法のメリットとしては，解析的に得られた数式を使うので数値計算誤差が出にくいことが挙げられる． 逆にデメリットとしては，手計算すれば分かるが計算が非常に面倒くさいことが挙げられる．もう少しモデルが複雑だったら発狂しかねない．

以下は解析的に求めたヤコビアンによるニュートン法の描画． きちんとそれっぽい位置に収束していることが分かる．(収束の様子を丁寧に見たかったので新しい仮定解を$\boldsymbol{x}_0=\boldsymbol{x}_0+0.1\Delta \boldsymbol{x}$としている)

In [65]:
## analysis 
def makeJx(sx):
    def Jacobian(usrx):
        x,y=usrx[0,0],usrx[1,0]
        r1,r2,r3=f1(usrx),f2(usrx),f3(usrx)
        j11=-(sx[0,0]-x)/r1
        j12=-(sx[1,0]-y)/r1
        j21=-(sx[0,1]-x)/r2
        j22=-(sx[1,1]-y)/r2
        j31=-(sx[0,2]-x)/r3
        j32=-(sx[1,2]-y)/r3
        return np.array([[j11,j12],[j21,j22],[j31,j32]])
    return Jacobian
#
J=makeJx(sx)
P = lambda A : np.linalg.solve( (A.T@A),A.T )
In [66]:
maxiter=50
x=array([[0],[0]])
xlist=np.zeros((2,maxiter))
for i in np.arange(0,maxiter):
    xlist[:,i]=x.ravel()
    dY=np.array([[r1-f1(x)],[r2-f2(x)],[r3-f3(x)]])
    dX=P(J(x))@dY
    x=x+dX*0.1

figure()
plot(sx[0,:],sx[1,:],'ob')
plot(xlist[0,:],xlist[1,:],'mo-')
plot_implicit(drawModel(sx1,r1),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx2,r2),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx3,r3),xrange=xrange,yrange=yrange)
grid()
xlim(xrange)
ylim(yrange)
Out[66]:
(-1, 4)

2.4数値微分でヤコビアンを求めて解く

仮定解$\boldsymbol{x}$が得られる度に，以下の式を計算して$\boldsymbol{x}_0$近傍のヤコビアンを得る方法である．

$$ \boldsymbol{J}_0 = \left( \begin{array}{c} \frac{f_1(\boldsymbol{x}+\boldsymbol{h}_x)-f_1(\boldsymbol{x})}{\boldsymbol{h}_x}&\frac{f_1(\boldsymbol{x}+\boldsymbol{h}_y)-f_1(\boldsymbol{x})}{\boldsymbol{h}_y}\\ \frac{f_2(\boldsymbol{x}+\boldsymbol{h}_x)-f_2(\boldsymbol{x})}{\boldsymbol{h}_x}&\frac{f_2(\boldsymbol{x}+\boldsymbol{h}_y)-f_2(\boldsymbol{x})}{\boldsymbol{h}_y}\\ \frac{f_3(\boldsymbol{x}+\boldsymbol{h}_x)-f_3(\boldsymbol{x})}{\boldsymbol{h}_x}&\frac{f_3(\boldsymbol{x}+\boldsymbol{h}_y)-f_3(\boldsymbol{x})}{\boldsymbol{h}_y}\\ \end{array} \right)$$

ただし，$\boldsymbol{h}_x=[e,0]^\mathrm{T},\boldsymbol{h}_y=[0,e]^\mathrm{T}$であり，$e$は適当な小さい数である．

In [67]:
### numerical yacob
def fprime2(fn):
    def fp(x):
        accr=0.001
        fpx=(fn(x+array([[accr],[0]]))-fn(x))/accr
        fpy=(fn(x+array([[0],[accr]]))-fn(x))/accr
        return array([[fpx,fpy]])
    return fp
fp1=fprime2(f1)
fp2=fprime2(f2)
fp3=fprime2(f3)
In [68]:
x=array([[0],[0]])
xlist=np.zeros((2,maxiter))
for i in np.arange(0,maxiter):
    xlist[:,i]=x.ravel()
    NJac=np.r_[fp1(x),fp2(x),fp3(x)]
    dY=np.array([[r1-f1(x)],[r2-f2(x)],[r3-f3(x)]])
    dX=P(NJac)@dY
    x=x+dX*0.1


figure()
plot(sx[0,:],sx[1,:],'ob')
plot(xlist[0,:],xlist[1,:],'mo-')
plot_implicit(drawModel(sx1,r1),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx2,r2),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx3,r3),xrange=xrange,yrange=yrange)
grid()
xlim(xrange)
ylim(yrange)
Out[68]:
(-1, 4)

この方法のメリットは面倒くさい偏微分計算をしなくていいこと． デメリットは，ヤコビアンがどんな数式か分からないので，計算が不安定になる条件とかが不明なこと(要出典)

おまけ

あまりにも解から遠い初期値で始めると局所解に陥る.

その境目はだいたい衛星同士を結んだ線．この線よりも上側を初期値にすると衛星よりも上側にある解に収束してしまう．

In [69]:
xrange=[-2,5]
yrange=[-2,5]
xq1,xq2=np.meshgrid(linspace(xrange[0],xrange[1],20),linspace(yrange[0],yrange[1],40))
xq=np.c_[xq1.ravel(),xq2.ravel()]

dxq=[]
for col in xq:
    x=vec(col)
    dY=np.array([[r1-f1(x)],[r2-f2(x)],[r3-f3(x)]])
    dX=P(J(x))@dY 
    dxq.append(dX)
dxq=array(dxq).squeeze(2)


figure()
plot(sx[0,:],sx[1,:],'ob')
plot_implicit(drawModel(sx1,r1),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx2,r2),xrange=xrange,yrange=yrange)
plot_implicit(drawModel(sx3,r3),xrange=xrange,yrange=yrange)
xq,dxq=xq.T,dxq.T
quiver(xq[0,:],xq[1,:],dxq[0,:],dxq[1,:]\
,angles='xy',color='m')
grid()
axis('equal')
xlim(xrange)
ylim(yrange)
Out[69]:
(-2, 5)
In [ ]: